Задача этого проекта — научиться предсказывать количество поездок в ближайшие часы в каждом районе Нью-Йорка. Для того, чтобы её решить, сырые данные необходимо агрегировать по часам и районам. Агрегированные данные будут представлять собой почасовые временные ряды с количествами поездок из каждого района. Похожие задачи возникают на практике, если вам необходимо спрогнозировать продажи большого количества товаров в большом количестве магазинов, объём снятия денег в сети банкоматов, посещаемость разных страниц сайта и т.д.
Авторы проекта предложили использовать инструмент на выбор. Т.к. для меня это второй проект, то я решил его сделать на R. Ввиду того, что в специализации используется Python, я постарался снабдить код подробными комментариями, чтобы (при желании) было легче разобраться.
Для работы потребуются следующие библиотеки:
library(leaflet) # интерактивные карты
library(scales) # в помощь графикам
library(ggplot2) # графики
library(ggmap) # геоданные с Гугла
library(magrittr) # c'est ne pas une pipe
library(data.table) # работа с таблицами данных
# Sys.setlocale('LC_ALL','utf-8') # если неверно отображается кириллица
Информация о версиях библиотек и системе:
sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] data.table_1.10.1 magrittr_1.5 ggmap_2.7
## [4] ggplot2_2.2.0.9000 scales_0.4.1 leaflet_1.0.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.8 knitr_1.15.1 maps_3.1.1
## [4] munsell_0.4.3 lattice_0.20-34 geosphere_1.5-5
## [7] colorspace_1.3-1 R6_2.2.0 rjson_0.2.15
## [10] jpeg_0.1-8 dplyr_0.5.0 stringr_1.1.0
## [13] plyr_1.8.4 tools_3.3.2 grid_3.3.2
## [16] gtable_0.2.0 png_0.1-7 DBI_0.5-1
## [19] htmltools_0.3.5 yaml_2.1.14 lazyeval_0.2.0
## [22] rprojroot_1.1 digest_0.6.10 assertthat_0.1
## [25] tibble_1.2 reshape2_1.4.2 mapproj_1.2-4
## [28] bitops_1.0-6 htmlwidgets_0.8 evaluate_0.10
## [31] rmarkdown_1.2 sp_1.2-3 stringi_1.1.2
## [34] RgoogleMaps_1.4.1 backports_1.0.4 proto_1.0.0
Просуммируйте общее количество поездок такси из каждой географической зоны и посчитайте количество ячеек, из которых в мае не было совершено ни одной поездки.
# загрузка
regions <- fread('data_in/regions.csv') # данные о регионах
dt_agg <- fread('data_out/2016_05.csv')
# суммирование по регионам
dt_month_regions <- dt_agg[, .(n = sum(n)), by = .(region)]
dt_month_regions[n == 0, .N]
## [1] 1283
Таким образом из 1283 регионов не было совершено ни одной поездки.
##########################
# Константа с геоданными #
##########################
# Нью-Йорк вписан в прямоугольник от -74.25559 до -73.70001 градусов долготы и от 40.49612 до 40.91553 широты
NY <- list(lon = c(-74.25559, -73.70001), lat = c(40.49612, 40.91553))
NY$CELLS <- 50 # число ячеек по широте или долготе
NY$lon_step <- (NY$lon[2] - NY$lon[1]) / NY$CELLS # шаг по долготе
NY$lat_step <- (NY$lat[2] - NY$lat[1]) / NY$CELLS # шаг по широте
NY$lon_breaks <- seq(NY$lon[1], NY$lon[2], by = NY$lon_step) # границы регионов по долготе
NY$lat_breaks <- seq(NY$lat[1], NY$lat[2], by = NY$lat_step) # границы регионов по широте
Поставьте на карте точку там, где находится Эмпайр-Стейт-Билдинг.
Для отображения статических карт будем пользоваться пакетами ggmap и ggplot2.
# координаты Empire State Building
# geocode('Empire State Building')
# введем координаты вручную
ESB <- list(lon = -73.98566, lat = 40.74844)
# скачиваем карту с маркером Эмпайр-Стейт-Билдинг
NY$map <- get_googlemap(c((NY$lon[1] + NY$lon[2]) / 2,
(NY$lat[1] + NY$lat[2]) / 2),
zoom = 10, scale = 2,
size = c(640, 640),
markers = data.frame(lon = ESB$lon, lat = ESB$lat)
)
# отображаем и обрезаем карту по границам регионов
ggmap(NY$map, extent = 'device') +
scale_x_continuous(limits = c(NY$lon[1], NY$lon[2])) +
scale_y_continuous(limits = c(NY$lat[1], NY$lat[2])) +
labs(title = 'Карта Нью-Йорка',
subtitle = 'Эмпайр-Стейт-Билдинг отмечен маркером')
Поверх статической карты Нью-Йорка визуализируйте данные о поездках из каждой ячейки так, чтобы цветовая шкала, в которую вы окрашиваете каждую ячейку, показывала суммарное количество поездок такси из неё.
Для начала посмотрим на распределение общего числа поездок по регионам за месяц. Для наглядности отобразим квадратный корень и логарифм этой величины:
# гистограмма n
dt_month_regions[n>0, .(region, n, sqrt_n = sqrt(n), log_n = log(n+1))] %>%
melt(measure.vars = 2:4, id.vars = 1) %>%
ggplot(aes(value, ..ncount..)) +
geom_histogram(bins = 15) +
facet_wrap(~variable, scales = 'free_x') +
labs(title = 'Гистограммы числа поездок',
subtitle = 'n, sqrt(n) и log(n+1)',
x = 'Значение', y = 'Норм. частота') +
scale_y_continuous(breaks = c(0, 0.5, 1))
Очевидно, что если не преобразовывать число поездок, то практически вся карта будет раскрашена одноцветно, т.к. много нулей и маленьких значений. Поэтому для раскраски каждого региона воспользуемся логарифмическим преобразованием, а чтобы ячейки с небольшими значениями имели меньший цветовой вес, пусть прозрачность определяется корнем величины. Регионы с нулевым значением отрисовывать не будем. В этом случае, тепловая карта суммарного числав поездок будет выглядеть так:
# добавляем каждому региону его границы
dt_month_regions <- cbind(dt_month_regions, regions[, 2:5])
# скачиваем чистую карту
NY$map <- get_googlemap(c((NY$lon[1] + NY$lon[2]) / 2,
(NY$lat[1] + NY$lat[2]) / 2),
zoom = 10, scale = 2,
size = c(640, 640)
)
# отображаем карту
ggmap(NY$map, extent = 'device') +
# рисуем квадратики, исключаем нулевые
geom_rect(data = dt_month_regions[n > 0],
aes(xmin = west, xmax = east, ymin = south, ymax = north,
fill = log(n), alpha = sqrt(n)),
size = 0.1, color = 'gray60', inherit.aes = FALSE) +
# обрезаем границы
scale_x_continuous(limits = c(NY$lon[1], NY$lon[2])) +
scale_y_continuous(limits = c(NY$lat[1], NY$lat[2])) +
# заливка квадратиков
scale_fill_distiller(palette = "OrRd", breaks = pretty_breaks(n = 6),
direction = 1) +
# прозрачность квадратиков в зависимости от n
scale_alpha(range=c(0.2, 0.5)) +
# убираем легенды
theme(legend.position="none") +
# заголовок
labs(title = 'Тепловая карта поездок за месяц')
На карту попали регионы, откуда поездок быть не может (бухты, река).
Вставьте интерактивную карту Нью-Йорка — такую, которую можно прокручивать и увеличивать. Поставьте метку там, где находится статуя свободы.
Для интерактивных карт воспользуемся пакетом leaflet.
# координаты Empire State Building
# geocode('Statue of Liberty')
# введем координаты вручную
SoL <- data.frame(lon = -74.0445, lat = 40.68925, name = 'Statue of Liberty')
# отображаем карту
leaflet(SoL, width = '100%') %>% addTiles() %>%
addMarkers(~lon, ~lat, popup = ~name) %>%
setView(SoL$lon, SoL$lat, zoom = 12)
Нарисуйте на интерактивной карте Нью-Йорка ячейки так, чтобы их цвет показывал среднее за месяц количество поездок такси в час из этой зоны.
Прежде всего необходимо рассчитать среднее количество поездок в час для каждого региона. Попробуем отобразить “горячую карту” так же, как и статичную карту.
Рассчитаем прозрачность как корень из числа поездок и отмасштабируем множество этих значений на отрезок [0.2; 0.6].
Для получения цвета для каждой зоны, создадим палитру из 50 оттенков (как в предыдущей тепловой карте). Среднее число поездок из каждого региона прологарифмируем, а все множество логарифмов разобьем на 50 интервалов, сопоставив таким образом каждому логарифму среднего числа поездок свой оттенок из палитры. Для удобства добавим текст по клику для квадратиков, а для самой карты предусмотрим возможность отключения как слоя. По аналогии с предыдущей картой, не будем отображать регионы, среднее значение которых меньше 0.1.
# рассчитываем среднее количество поездок
dt_agg_avr <- dt_agg[, .(avr_n = mean(n)), by = .(region)]
# добавляем границы квадратиков
dt_agg_avr <- cbind(dt_agg_avr, regions[,2:5])
# палитра для квадратиков
shades <- 50 # 50 оттенков красно-оранжевого
pal_OrRd <- seq(0, 1, length.out = 50) %>%
(brewer_pal(palette = 'OrRd', direction = -1)(9)[1:5] %>%
rev %>% gradient_n_pal)
# прологарифмируем средние значения и разобьем их на 50 градаций
# сопоставляем цвет из палитры каждому полученному числу
dt_agg_avr <- dt_agg_avr[,
fill_color := pal_OrRd[log(avr_n + 1) %>% cut(breaks = shades) %>%
as.numeric()]]
# добавляем признак прозрачности, ~ sqrt(n)
dt_agg_avr <- dt_agg_avr[, alpha := sqrt(avr_n + 1) %>% rescale(c(0.2, 0.6))]
# отрисуем карту
leaflet(dt_agg_avr[avr_n > 0.1], width = '100%') %>% addTiles() %>%
# квадратики
addRectangles(lng1 = ~west, lng2 = ~east, lat = ~south, lat2 = ~north,
# границы квадратиков
weight = 1, color = 'gray', opacity = 4/10,
# заливка квадратиков
fillOpacity = ~alpha, fillColor = ~fill_color,
# описание квадратика
popup = ~paste('region:', region, 'avr_n:', round(avr_n, 1)),
# для управления слоями
group = 'heatmap') %>%
# чтобы квадратики можно было убирать
addLayersControl(
overlayGroups = c("heatmap"),
options = layersControlOptions(collapsed = FALSE))
Чтобы не выбирать из всех 2500 ячеек вручную, отфильтруйте ячейки, из которых в мае совершается в среднем меньше 5 поездок в час. Посчитайте количество оставшихся. Проверьте на карте, что среди этих ячеек нет таких, из которых поездки на самом деле невозможны.
После того, как уберем зоны, откуда было в среднем было меньше 5 поездок в час, останется 102 региона:
dt_agg_avr[avr_n >= 5, .N]
## [1] 102
Отобразим и выделим их на карте. Включая и отключая слой с регионами можно убедиться, что из всех этих зон возможны поездки.
# отрисуем карту
leaflet(dt_agg_avr[avr_n >= 5], width = '100%') %>% addTiles() %>%
# квадратики
addRectangles(lng1 = ~west, lng2 = ~east, lat = ~south, lat2 = ~north,
# границы квадратиков
weight = 1, color = 'gray', opacity = 3/4,
# заливка квадратиков
fillOpacity = 1/2, fillColor = 'maroon',
# описание квадратика
popup = ~paste('region:', region, 'avr_n:', round(avr_n, 1)),
# для управления слоями
group = 'heatmap') %>%
# чтобы квадратики можно было убирать
addLayersControl(
overlayGroups = c("heatmap"),
options = layersControlOptions(collapsed = FALSE))
Ноутбук опубликован на rstudioconnect.com. Файлы проекта также можно найти и на гитхабе.